# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# File Name          : run_paper_analysis_noPS.R
# Programmer Name    : Luis F. Campos
#                     luiscampos@g.harvard.edu
#
# Purpose            : This file contains code to run the analysis presented 
#					   in Section 7 of "Worth Weighting? How to Think 
#					   About and Use Sample Weights in Survey Experiments"
#					   In particular it analyzes the survey experiments 
#					   with no post-stratification
#
# Input              : 'survey_exp_toolkit.R', 
#					   'function_analysis_multilevel_treatment.R'
#					   'run_paper_data_analysis_sub.R'
#					   'Survey_Experiment_Data.Rda'
#
# Output             : final_results_no_PS.csv
# 
# References         : None
#
#
# Platform           : R
# Version            : v3.3.0
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Replication Information:
#   We do expect differences when running this code on different systems 
#   (and with different seeds), but the differences all tend to be quite 
#   small and the general trends remain. If you would like to replicate 
#   the exact numbers and figures used in the article, in the readme you 
#   will find the session info. Please install and load all package 
#   versions (including the R version) to ensure the exact replication. 
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #



set.seed(2)
# Set current working directory
# setwd('')

# load basic functions
source("./_survey_exp_toolkit.R")
# load analysis function
source("./_function_analysis.R")
# load experiments for analysis 
source("./_paper_analysis_labels.R")


out_dir = './'


library( plyr )
library( xtable )
library( reshape2 )
library(foreign, pos=4)
library(Matching)
library(car)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Number of Bootstrat samples for SE Estimates
#  Change to smaller number to speed things up
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

B = 1000

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# Load Survey Experiment Data
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #


load(file = 'Survey_Experiment_Data.Rda')




# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 	Virginia
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

results_Virginia = list()
# remove PS label replicates and datasets (i.e. keep every other experiment label)
TMT_Virginia_noPS = TMT_Virginia[seq(1, nrow(TMT_Virginia), 2),]
dat_Virginia_noPS = dat_Virginia[seq(1, nrow(TMT_Virginia), 2)]

for(i in 1:nrow(TMT_Virginia_noPS)){
	r = TMT_Virginia_noPS[i,]
	print(r)
	dat = dat_Virginia_noPS[[i]]
	dat$noPS = '0'
	dat$noPS[which(as.character(dat[,r[6]]) %in% c('Republican', 'Democrat'))] = '1'
	results_Virginia[[i]] = do.analysis(dat = dat, w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = 'noPS', sub_level = '1', B = B, return.all = TRUE)
}


save(results_Virginia, file = paste(out_dir, 'TMT_Virginia_noPS.Rda', sep = ''))





# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 	CCES12_AGT_Matched_working
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #



results_CCES12 = list()
# remove PS label replicates (i.e. keep every other experiment label)
TMT_CCES12_noPS = TMT_CCES12[seq(1, nrow(TMT_CCES12), 2),]
dat_CCES12_noPS = dat_CCES12[seq(1, nrow(TMT_CCES12), 2)]

for(i in 1:nrow(TMT_CCES12_noPS)){
	r = TMT_CCES12_noPS[i,]
	print(r)
	dat = dat_CCES12_noPS[[i]]
	dat$noPS = '0'
	dat$noPS[which(as.character(dat[,r[6]]) %in% c('Republican', 'Democrat'))] = '1'
	results_CCES12[[i]] = do.analysis(dat = dat, w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = 'noPS', sub_level = '1', B = B, return.all = TRUE)
}

save(results_CCES12, file = paste(out_dir, 'TMT_CCES12_noPS.Rda', sep = ''))




# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# YouGov
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #


results_YouGov = list()
# remove PS label replicates (i.e. keep every other experiment label)
TMT_YouGov_noPS = TMT_YouGov[seq(1, nrow(TMT_YouGov), 2),]
dat_YouGov_noPS = dat_YouGov[seq(1, nrow(TMT_YouGov), 2)]


for(i in 1:nrow(TMT_YouGov_noPS)){
	r = TMT_YouGov_noPS[i,]
	print(r)
	dat = dat_YouGov_noPS[[i]]
	dat$noPS = '0'
	dat$noPS[which(as.character(dat[,r[6]]) %in% c('Republican', 'Democrat'))] = '1'
	results_YouGov[[i]] = do.analysis(dat = dat, w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = 'noPS', sub_level = '1', B = B, return.all = TRUE)
}

save(results_YouGov, file = paste(out_dir, 'TMT_YouGov_noPS.Rda', sep = ''))



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 2010_CCES_UCB
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

results_CCES10_UCB = list()
# remove PS label replicates (i.e. keep every other experiment label)
TMT_CCES10_UCB_noPS = TMT_CCES10_UCB[seq(1, nrow(TMT_CCES10_UCB), 2),]
dat_CCES10_UCB_noPS = dat_CCES10_UCB[seq(1, nrow(TMT_CCES10_UCB), 2)]

for(i in 1){
	r = TMT_CCES10_UCB_noPS
	print(r)
	dat = dat_CCES10_UCB_noPS[[i]]
	dat$noPS = '0'
	dat$noPS[which(as.character(dat[,r[6]]) %in% c('Republican', 'Democrat'))] = '1'
	results_CCES10_UCB[[i]] = do.analysis(dat = dat, w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = 'noPS', sub_level = '1', B = B, return.all = TRUE)
}

save(results_CCES10_UCB, file = paste(out_dir, 'TMT_CCES10_UCB_noPS.Rda', sep = ''))





# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# CCES10
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

results_CCES10_EIT = list()
# remove PS label replicates (i.e. keep every other experiment label)
TMT_CCES10_EIT_noPS = TMT_CCES10_EIT[seq(1, nrow(TMT_CCES10_EIT), 2),]
dat_CCES10_EIT_noPS = dat_CCES10_EIT[seq(1, nrow(TMT_CCES10_EIT), 2)]


for(i in 1:nrow(TMT_CCES10_EIT_noPS)){
	r = TMT_CCES10_EIT_noPS[i,]
	print(r)
	dat = dat_CCES10_EIT_noPS[[i]]
	dat$noPS = '0'
	dat$noPS[which(as.character(dat[,r[6]]) %in% c('Republican', 'Democrat'))] = '1'
	results_CCES10_EIT[[i]] = do.analysis(dat = dat, w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = 'noPS', sub_level = '1', B = B, return.all = TRUE)
}

save(results_CCES10_EIT, file = paste(out_dir, 'TMT_CCES10_EIT_noPS.Rda', sep = ''))



# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 	UCB_Follow
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

results_UCB_Follow = list()
# remove PS label replicates (i.e. keep every other experiment label)
TMT_UCB_Follow_noPS = TMT_UCB_Follow[seq(1, nrow(TMT_UCB_Follow), 2),]
dat_UCB_Follow_noPS = dat_UCB_Follow[seq(1, nrow(TMT_UCB_Follow), 2)]


for(i in 1:nrow(TMT_UCB_Follow_noPS)){
	r = TMT_UCB_Follow_noPS[i,]
	print(r)
	dat = dat_UCB_Follow_noPS[[i]]
	dat$noPS = '0'
	dat$noPS[which(as.character(dat[,r[6]]) %in% c('Republican', 'Democrat'))] = '1'
	results_UCB_Follow[[i]] = do.analysis(dat = dat, w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = 'noPS', sub_level = '1', B = B, return.all = TRUE)
}

save(results_UCB_Follow, file = paste(out_dir, 'TMT_UCB_Follow_noPS.Rda', sep = ''))




# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 2014_CCES_Matched: DO
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #


results_2014_CCES_matched = list()
# remove PS label replicates (i.e. keep every other experiment label)
TMT_2014_CCES_matched_noPS = TMT_2014_CCES_matched[seq(1, nrow(TMT_2014_CCES_matched), 2),]
dat_2014_CCES_matched_noPS = dat_2014_CCES_matched[seq(1, nrow(TMT_2014_CCES_matched), 2)]


for(i in 1:nrow(TMT_2014_CCES_matched_noPS)){
	r = TMT_2014_CCES_matched_noPS[i,]
	print(r)
	dat = dat_2014_CCES_matched_noPS[[i]]
	dat$noPS = '0'
	dat$noPS[which(as.character(dat[,r[6]]) %in% c('Republican', 'Democrat'))] = '1'

	results_2014_CCES_matched[[i]] = do.analysis(dat = dat, w = r[1], treatment = r[2], t0 = r[3], t1 = r[4], outcome = r[5], sub_var = 'noPS', sub_level = '1', B = B, return.all = TRUE)
}

save(results_2014_CCES_matched, file = paste(out_dir, 'TMT_2014_CCES_matched_noPS.Rda', sep = ''))





# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# 	Gather results together
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
out_dir = './'


load(file = paste(out_dir, 'TMT_Virginia_noPS.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_CCES12_noPS.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_YouGov_noPS.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_CCES10_UCB_noPS.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_2014_CCES_matched_noPS.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_CCES10_EIT_noPS.Rda', sep = ''))
load(file = paste(out_dir, 'TMT_UCB_Follow_noPS.Rda', sep = ''))


results = c(results_Virginia, results_CCES12, results_YouGov, results_CCES10_UCB, results_2014_CCES_matched, results_CCES10_EIT, results_UCB_Follow)

nms = rep(c('Virginia', 'CCES12', 'YouGov', 'CCES10_UCB', 'CCES2014', 'CCES10_EIT', 'UCB_Follow'), c(length(results_Virginia), length(results_CCES12), length(results_YouGov), length(results_CCES10_UCB), length(results_2014_CCES_matched), length(results_CCES10_EIT), length(results_UCB_Follow)))




res_SD = sapply(results, function(x){
	c(x$res[4,3]/x$res[7,3], x$dats['n'])
})


r.keep = c(1, 4, 7, 8)
morph = function(tab){
	nm = do.call('c', lapply(tab[r.keep,1], function(x) paste(x, c('.est', '.SE'), sep = '')))
	est = c(as.matrix(t(tab[r.keep,c(2, 3)])))
	names(est) = nm
	est
}

dats1 = lapply(results, function(x)x$dats)
res1 = lapply(results, function(x)morph(x$res))

res = data.frame(do.call('rbind', dats1), do.call('rbind', res1), survey = nms)




# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# CIs and Significance
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
res$tau.sate.L = res[,"tau.sate.est"] - 1.96*res[,"tau.sate.SE"]
res$tau.sate.U = res[,"tau.sate.est"] + 1.96*res[,"tau.sate.SE"]
res$tau.hh.L = res[,"tau.hh.est"] - 1.96*res[,"tau.hh.SE"]
res$tau.hh.U = res[,"tau.hh.est"] + 1.96*res[,"tau.hh.SE"]
res$tau.ps.hh.L = res[,"tau.ps.hh.est"] - 1.96*res[,"tau.ps.hh.SE"]
res$tau.ps.hh.U = res[,"tau.ps.hh.est"] + 1.96*res[,"tau.ps.hh.SE"]

res$tau.sate.sig = ((res$tau.sate.L > 0) | (res$tau.sate.U < 0))
res$tau.hh.sig = ((res$tau.hh.L > 0) | (res$tau.hh.U < 0))
res$tau.ps.hh.sig = ((res$tau.ps.hh.L > 0) | (res$tau.ps.hh.U < 0))

res$delta = (res[,"tau.sate.est"]-res[,"tau.hh.est"])/res[,"tau.star1.SE"]


write.csv(res, file = paste(out_dir, 'final_results_no_PS.csv', sep = ''))






